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Abstract — This paper offers a natural stochastic semantics of 
Networks of Priced Timed Automata (NPTA) based on races 
between components. The semantics provides the basis for 
satisfaction of Probabilistic Weighted CTL properties (PWCTL), 
conservatively extending the classical satisfaction of timed auto- 
mata with respect to TCTL. In particular the extension allows 
for hard real-time properties of timed automata expressible in 
TCTL to be refined by performance properties, e.g. in terms of 
probabilistic guarantees of time- and cost-bounded properties. A 
second contribution of the paper is the application of Statistical 
Model Checking (SMC) to efficiently estimate the correctness 
of non-nested PWCTL model checking problems with a desired 
level of confidence, based on a number of independent runs of 
the NPTA. In addition to applying classical SMC algorithms, 
we also offer an extension that allows to efficiently compare 
performance properties of NPTAs in a parametric setting. The 
third contribution is an efficient tool implementation of our result 
and applications to several case studies. 

I. Introduction 

Model Checking (MC) JTJ is a widely recognised approach 
to guarantee the correctness of a system by checking that any 
of its behaviors is a model for a given property. There are 
several variants and extensions of MC aiming at handling 
real-time and hybrid systems with quantitative constraints 
on time, energy or more general continuous aspects ||2]-||5). 
Within the field of embedded systems these formalisms and 
their supporting tools [6|-[9) are now successfully applied 
to time- and energy-optimal scheduling, WCET analysis and 
schedulability analysis. 

Compared with traditional approaches, a strong point of 
real-time model checking is that it (in principle) only requires 
a model to be applicable, thus extensions to multi-processor 
setting is easy. A weak point of model checking is the 
notorious problem of state-space explosion, i.e. the exponential 
growth in the analysis effort measured in the number of model- 
components. Another limitation of real-time model checking 
is that it merely provides - admittedly most important - 
hard quantitative guarantees, e.g. the worst case response time 
of a recurrent task under a certain scheduling principle, the 
worst case execution time of a piece of code running on a 
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particular execution platform, or the worst case time before 
consensus is reached by a real-time network protocol. In 
addition to these hard guarantees, it would be desirable in 
several situations to obtain refined performance information 
concerning likely or expected behaviors in terms of timing 
and resource consumption. In particular, this would allow to 
distinguish and select between systems that perform identically 
from a worst-case perspective. 

To illustrate our point consider the network of two priced 
timed automata in Fig.[T]modeling a competition between Axel 
and Alex both having to hammer three nails down. As can be 
seen by the representing Work-locations the time (-interval) 
and rate of energy-consumption required for hammering a nail 
depends on the player and the nail-number. As expected Axel 
is initially quite fast and uses a lot of energy but becomes 
slow towards the last nail, somewhat in contrast to Alex. To 
make it an interesting competition, there is only one hammer 
illustrated by repeated competitions between the two players 
in the Ready-locations, where the slowest player has to 
wait in the Idle-location until the faster player has finished 
hammering the next nail. Interestingly, despite the somewhat 
different strategy applied, the best- and worst-case completion 
times are identical for Axel and Alex: 59 seconds and 150 
seconds. So, there is no difference between the two players and 
their strategy, or is there? Assume that a third person wants to 
bet on who is the more likely winner - Axel or Alex - given 
a refined semantics, where the time-delay before performing 
an output is chosen stochastically (e.g. by drawing from a 
uniform distribution). Under such a refined semantics there is 
a significant difference between the two players. In Fig. [2^) the 
probability distributions for either of the two players winning 
before a certain time is given. Though it is clear that Axel 
has a higher probability of winning than Alex (59% versus 
41%), however declaring the competition a draw if it has not 
finished before 50 seconds actually makes Alex the more likely 
winner. Similarly, Fig. |2j>) illustrates the probability of either 
of the two players winning given an upper bound on energy. 
With an unlimited amount of energy, clearly Axel is the most 
likely winner, whereas limiting the consumption of energy to 
maximum 52 "energy-units" gives Alex an advantage. 

As a first contribution of this paper we propose a stochastic 
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Figure 1: 3-Nail Hammering Game between Axel and Alex. 




semantics for Priced Timed Automata (PTA), whose clocks 
can evolve with different rates, whil^U being used with no re- 
strictions in guards and invariants. Networks of PTAs (NPTA) 
are created by composing PTAs via input and output actions. 
The model is as expressive as linear hybrid automata 0, 
making even the reachability problem undecidable. More pre- 
cisely, we define a natural stochastic semantics for networks of 
NPTAs based on races between components. We shall observe 
that such race can generate arbitrarily complex stochastic 
behaviors from simple assumptions on individual components. 
While fully stochastic semantics have already been proposed 
for timed systems iflOl . ifTTl . we are the first to consider 
networks of timed and hybrid systems. Other related work 
includes the very rich framework of stochastic timed systems 
of MoDeST fl2l . Here, however, general hybrid variables 
are not considered and parallel composition does not yield 
fully stochastic models. For the notion of probabilistic hybrid 
systems considered in lfl3ll the choice of time is resolved 
non-deterministic ally rather than stochastically as in our case. 
Moreover, based on the stochastic semantics, we are able 
to express refined performance properties, e.g. in terms of 
probabilistic guarantees of time- and cost-bounded properties. 

To allow for the efficient analysis of probabilistic perform- 
ance properties - despite the general undecidability of these - 

'in contrast to the usual restriction of priced timed automata [4|, |5| 



we propose to work with Statistical Model Checking (SMC) 
[14 1, [ 15 1, an approach that has recently been proposed as an 
alternative to avoid an exhaustive exploration of the state-space 
of the model. The core idea of the approach is to monitor 
some simulations of the system, and then use results from 
the statistic area (including sequential hypothesis testing or 
Monte Carlo simulation) in order to decide whether the system 
satisfies the property or not with some degree of confidence. 
By nature, SMC is a compromise between testing and classical 
model checking techniques. 

Thus, as a second contribution, we provide an efficient 
implementation of existing SMC algorithms that we use for 
checking the correctness of NPTAs with respect to cost- 
constrained temporal logic. The series of algorithms we im- 
plement includes a version of the sequential hypothesis test 
by Wald [16| as well as a quantitative approach [17|. Our 
implementation relies on a new efficient algorithm for genera- 
ting runs of NPTAs in a random manner. In addition, we also 
propose another SMC algorithm to compare the performances 
of two properties without computing their probability. This 
problem, which is far beyond the scope of existing time model 
checking approaches, can be approximated with an extension 
of the sequential hypothesis testing. In addition to be the first 
to apply such extension in the context of formal verification, 
we also propose a new variant that allows to reuse existing 
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results in parallel when comparing the properties on different 
timed bounds. 

Finally, one of the most interesting contribution of our 
work takes the form of a series of new case studies that 
are analyzed with a new stochastic extension of Uppaal 
|18|. Particularly, we show how our approach can be used 
to resolve scheduling problems. Such problems are defined 
using Duration Probabilistic Automata (DPA) lfl9l . a new and 
natural model for specifying list of tasks and shared resources. 
We observe that our approach is not only more general, but 
also much faster than the hypothesis testing engine recently 
implemented in the Prism toolset. Our work thus presents 
significant advances in both the modeling and the efficient 
verification of network of complex systems. 

Related work. Some works on probabilistic semantics 
of timed automata have already been discussed above. 
Simulation-based approaches such as Monte Carlo have been 
in use since decades, however the use of simulation and hypo- 
thesis testing to reason on formal models is a more recent ad- 
vance. First attempts to apply hypothesis testing on stochastic 
extension of Hennessy-Milner logic can be found in [20]. 
In fl4l . ||2T1 . Younes was the first to apply hypothesis testing 
to stochastic systems whose properties are specified with 
(bounded) temporal logic. His approach is implemented in the 
Ymer toolset |22] and can be applied on time-homogeneous 
generalized semi-Markov processes, while our semantics ad- 
dresses the composition of stochastic systems allowing to 
compose a global system from components and reason about 
communication between independent processes. In addition to 
Younes work we explore continuous-time features, formalize 
and implement Wald's ideas where the probability comparison 
can be evaluated on NPTA processes. In a recent work l23l . 
Zuliani et al. extended the SMC approach to hybrid systems. 
Their work is a combination of ll24l and 11251 based on 
Simulink models (non-linear hybrid systems), whereas our 
method is specialised to networks of priced timed automata 
where model-checking techniques can be directly applicable 
using the same tool suite. In addition we provide means 
of comparing performances without considering individual 
probabilities. Finally, a very recent work [26] proposes par- 
tial order reduction techniques to resolve non-determinism 
between components rather than defining a unique stochastic 
distribution on their product behaviors. While this work is of 
clear interest, we point out that the application of partial order 
may considerably increase the computation time and for some 
models partial orders cannot resolve non-determinism, espe- 
cially when considering continuous time l27l . Other works on 
SMC can be found in |28l. l29l. 

II. Network of Priced Timed Automata 

We consider the notion of Networks of Priced Timed Auto- 
mata (NPTA), generalizing that of regular timed automata (TA) 
in that clocks may have different rates in different locations. 
In fact, the expressive power (up to timed bisimilarity) of 
NPTA equals that of general linear hybrid automata (LHA) 



0, rendering most problems - including that of reachability 
- undecidable. 

Let X be a finite set of variables, called c/ocfcQ. A clock 
valuation over X is a mapping v : X — > R>o, where R>o is 
the set of nonnegative reals. We write K.> for the set of clock 
valuations over X. Let r : X — > N be a rate vector, assigning 
to each clock of X a rate. Then, for v E R> and d E M>o 
a delay, we write v + r ■ d for the clock valuation defined 
by [y + r ■ d)(x) = v(x) + r[x) ■ d for any clock x E X. 
We denote by N x the set of all rate vectors. If Y C X, 
the valuation v\Y] is the valuation assigning when x E Y 
and v(x) when x ^ Y . An upper bounded (lower bound) 
guard over X is a finite conjunction of simple clock bounds 
of the form x ~ n where x E X, n E N, and ^E {<, <} 
(~G {>,>}) We denote by U{X) (C{X) the set of upper 
(lower) bound guards over X, and write v \= g whenever v is 
a clock valuation satisfying the guard g. Let E = Ej tfcl E D be 
a disjoint sets of input and output actions. 

Definition 1 A Priced Timed Automaton (PTA) is a tuple A = 
(L, £q, X, E, E, R, I) where: (i) L is a finite set of locations, 
(ii) £q E L is the initial location, (iii) X is a finite set of clocks, 
(iv) E = Ej W S is a finite set of actions partitioned into inputs 
(£,) and outputs (E D ), (v) E C L x £(X) x E x 2 X x L is a 
finite set of edges, (vi) R : L — > N x assigns a rate vector to 
each location, and (viii) I : L — > U(X) assigns an invariant 
to each location. 

The semantics of NPTAs is a timed labelled transition system 
whose states are pairs (I, v) G L x R> with v |= 1(1), and 

whose transitions are either delay (£, v) — — > (£, v') with d E 
R> Q and v' = v + R(£) ■ d, or discrete (£, u) (£', v') if 
there is an edge (£, g, a, Y, £') such that v |= g and v' = u[Y]. 
We write (£, v) ~> (£! ', v') if there is a finite sequence of delay 
and discrete transitions from (£,v) to (£' ,v'). 

a) Networks of Priced Timed Automata: Following 
the compositional specification theory for timed systems in 
ll30ll . we shall assume that NPTAs are: (l)[Input-enabled:] for 
all states (£, v) and input actions t E (£, v) — >, and 
(2) [Deterministic:] for all states (£, v) and actions a G E, 
whenever (£, v) (£', v') and (£, v) {£" , v") then £' = 
£" and v' = v". 

Whenever A> = ;/.'. .Y /v'. /<"./' : (J = 1 . . .n) are 
NPTA, they are composable into a closed network iff their 
clock sets are disjoint (X^ n X k = when j ^ k), they have 
the same action set (E = E J = E fc for all j, k), and their 
output action-sets provide a partition of E (££ n E* = for 
j ^ k, and E = Uj£j). For a E E we denote by c(a) the 
unique j with a E E- 7 . 

Definition 2 Let A j = ■/.'. .Y'.X. /•.'. //'./': (j = l...n) 
be composable NPTAs. Then the composition (Ai . . . | A n ) 
is the NPTA A = (L, X, E, E, R, L) where (i) L = XjL j , 

2 We will (mis)use the term "clock" from timed automata, though in the 
setting of NPTAs the variables in X are really general real-valued variables. 
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(ii) X = UjX?, (iii) R{£){x) = RP(£ j )(x) when x G X-?', (iv) 
/(£) = r\jl[£ j ), and (v) (£, Hjgj, a, Uj-r,, ■£') G £ whenever 
{tjiQj^^jJ'j) G i? 3 for j = 1 . . .71. 
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Figure 3: Four composable NPTAs: A, B and T; A, B r and 
T; and AB and T. 

Example 1. Let A, £?, T and AB be the priced timed automata 
depicted in Fig. |J[ Then A, B and T are composable as 
well as AB and T. In fact the composite systems (A\B\T) 
and (AB\T) are timed (and priced) bisimilar, both having the 
transition sequence: 

((Ao,B„,T o ),[x = 0,y = 0,C = 0]) AA 
((i4 1 ,B ,T 1 ),ta; = l,i/ = l,C = 4]) -^A 
((j4i,Bi,T 3 ),[a: = 2,y = 2,C = 6]) I 

demonstrating that the final location T3 of T is reachable with 
cost 6. 

III. Probabilistic Semantics of NPTA 

Continuing Example 1 we may realise that location T3 of 
the component T is reachable within cost to 6 and within 
total time and 2 in both (A\B\T) and (AB\T) depending on 
when (and in which order) A and B {AB) chooses to perform 
the output actions a! and b\. Assuming that the choice of these 
time-delays is governed by probability distributions, we will 
in this section define a probability measure over sets of infinite 
runs of networks of NPTAs. 

In contrast to the probabilistic semantics of timed automata 
in iflOll . ifTTl our semantics deals with networks and thus with 
races between components. Let A' = (L 3 , X 3 , S, E J , R? , P) 
(j = 1 ... 71) be a collection of composable NPTAs. Under 
the assumption of input-enabledness, disjointness of clock sets 
and output actions, states of the the composite NPTA A. = 
(A\ I ... I A n ) may be seen as tuples s = (s%, . . . , s n ) where 
Sj is a state of A 3 , i.e. of the form (£, v) where I G L 3 
and v G K>q- Our probabilistic semantics is based on the 
principle of independency between components. Repeatedly 
each component decides on its own - based on a given delay 
density function and output probability function - how much 
to delay before outputting and what output to broadcast at that 
moment. Obviously, in such a race between components the 
outcome will be determined by the component that has chosen 
to output after the minimum delay: the output is broadcast and 
all other components may consequently change state. 

3 it is assumed that all components are completed with looping input 
transitions, where these are missing. 



b) Probabilistic Semantics of NPTA Components: 

Let us first consider a component A and let St-* denote the 
corresponding set of states. For each state s = (£, v) of A 
we shall provide probability distributions for both delays and 
outputs. 

The delay density function p s over delays in M>o will be 
either a uniform or an exponential distribution depending on 
the invariant of £. Denote by Eg the disjunction of guards 
g such that (£, g, o, — ,— ) G E 3 for some output o. Denote 
by d(£, v) the infimum delay before enabling an output, i.e. 
d(£, v) = inf {d G K> : v + RP ■ d \= E e }, and denote by 
D(£, v) the supremum delay, i.e. D(£, v) = sup{d G M>o : 
v + W ■ d |= P{£)}. If D(£, v) < 00 then the delay density 
function fi s is a uniform distribution on [d(£, v), D(£, u)]. 
Otherwise - that is P(£) does not put an upper bound on 
the possible delays out of s - the delay density function 
/i s is an exponential distribution with a rate P(£), where 
P : L- ' — » R>o is an additional distribution rate component 
added to the NPTA A 3 . For every state s — (£, v), the output 
probability function ^ s over S 3 is the uniform distribution over 
the set {o : (£, g,o,—,—) G E 3 A v |= g} whenever this set is 
non-empty 0. We denote by s° the state after the output of o. 
Similarly, for every state s and any input action l, we denote 
by s c the state after having received the input l. 

c) Probabilistic Semantics of Networks of NPTA: 
We shall now see that while the stochastic semantics of each 
PTA is rather simple (but quite realistic), arbitrarily complex 
stochastic behavior can be obtained by their composition. 

Reconsider the closed network A = (Ai . . . | A n ) with a 
state space St = Sti x • • • x St n . For s = (si, . . . , s n ) G St and 
a\di . . . cik G S* we denote by 7r(s, a\a,i . . . a^) the set of all 
maximal runs from s with a prefix tiai^Q^ • • • ifcCtfc for some 
ti,...,t n G K>o, that is runs where the i'th action has 
been outputted by the component A c ^ a .y We now inductively 
define the following measure for such sets of runs: 

P^.(7r(s, aia 2 ■ ..a n )) = 

It>o /*•. W ' ( Il^c I T >t Vsj ( T ) dT ) ■ TV (ai) ' 
P A (vr(s*) a S a 2 ...a n )) dt 

where c = c(ai), and as base case we take P4.(7r(s),e) = 1. 

This definition requires a few words of explanation: at the 
outermost level we integrate over all possible initial delays 
t. For a given delay t, the outputting component c = c(ai) 
will choose to make the broadcast at time t with the stated 
density. Independently, the other components will choose to 
a delay amount, which - in order for c to be the winner - 
must be larger than t; hence the product of the probabilities 
that they each make such a choice. Having decided for making 
the broadcast at time t, the probability of actually outputting 
ai is included. Finally, in the global state resulting from all 
components having delayed t time-units and changed state 
according to the broadcasted action a\ the probability of runs 

4 otherwise a specific weight distribution can be specified and used instead. 
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(a) (b) 
Figure 4: Cumulative probabilities for time and cost-bounded 
reachability of T 3 . 

according to the remaining actions 0,2 . . . a n is taken into 
account. 

d) Logical Properties: Following [31 |, the measure 
may be extended in a standard and unique way to the a- 
algebra generated by the sets of runs (so-called cylinders) 
7r(s,aia2 . . . a n ). As we shall see this will allow us to give 
proper semantics to a range of probabilistic time- and cost- 
constrained temporal properties. Let A be a NPTA. Then we 
consider the following non-nested PWCTL properties: 



"(Oc<c^) ~p I P(nc< c <^) ~p 



where C is an observer clock (of A), if a state-property (wrt. 
.4) , {<, <, =, >, >}, and p e [0, 1]. For the semantics 
let A* be the modification of A, where the guard C < c has 
been conjoined to the invariant of all locations and an edge 
(£,(p,Otp,®,£) has been added to all edges £, where o v is a 
new output action. Then: 



A\=V(O c <cV) ~P iff Px'( (J 

7r(s , OO, 



crGS* 



■ p 



which is well-defined since the a-algebra on which P_4* is 
defined is closed under countable unions and finite intersec- 
tions. To complete the semantics, we note that P(nc<cV) 
is equivalent to (1 — p) ~ ¥{<>c<c^<-p)- 
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Example 1 Reconsider the Example of Fig. [3] Then it can be 
shown that (A\B\T) \= P(0 4 < 2 T 3 ) = 0.75 and (A\B\T) \= 
^(Oc<6T 3 ) = 0.75, whereas (AB\T) \= P(O t < 2 T 3 ) = 0.50 
and (AB\T) \= P(O c < 6 T 3 ) = 0.50. Fig. g] gives a time- 
and cost-bounded reachability probabilities for (A|_B|T) and 
(AB\T) for a range of bounds. Thus, though the two NPTAs 
satisfy the same WCTL properties, they are obviously quite 
different with respect to PWCTL. The NPTA B r of Fig. [3]is a 
variant of B, with the uniform delay distribution enforced by 
the invariant y < 2 being replaced by an exponential distribu- 
tion with rate \. Here (A\B r \T) satisfies P(O f < 2 T 3 ) k, 0.41 
and P(O c < 6 T 3 ) w 0.49. 

IV. Statistical Model Checking for NPTA 

As we pointed out, most of model checking problems for 
NPTAs and PWCTL (including reachability) are undecidable. 
Our solution is to use a technique that approximates the an- 
swer. We rely on Statistical Model Checking (SMC) 03], fI31 . 
that is a series of simulation-based techniques that generate 
runs of the systems, monitor them, and then use algorithms 
from statistics to get an estimate of the entire system. At the 



Algorithm 1: Random run for a NPTA-network A 

function RRa((£o, vq),C, c) 

1 run \— i£,v) := tail(ruri) := {£0,1^0) 

2 while v (C) < c do 

3 for i = 1 to \£\ do di := delay(fi^ e . iV .- ) ) 

4 d := min 1 <i<\ i \(di) 

5 if d = +00 V v(C) + d * R{£){C) > c then 

6 d:= (u(G) -c)/R(£)(G) 

return run® —>■ (£, v + d * R(£)) 

end 
else 

pick k such that d^ — d; Vd := v + d * R(£) 
pick 4 -^^> E' k with g{v d ) 



9 
10 
11 



run := runt 



end 

(/,"):= 
end 

return run 



$(t,v d )^(t[l' k /l k ],[r*0](v d )) 



tail(run) 



heart of any SMC approach, there is an algorithm used to 
generate runs of the system following a stochastic semantics. 
We propose such an algorithm for NPTAs corresponding to 
the stochastic semantics proposed in Section [TTIJ Then, we 
recap existing statistic algorithms, providing the basis for a 
first SMC algorithm for NPTAs. 

e) Generating Runs of NPTA: SMC is used for prop- 
erties that can be monitored on finite runs. Here, we propose 
an algorithm that given an NPTA generates a random run up 
to a cost bound c (with time bounds being a simple case) 
of an observer clock C. A run of a NPTA is a sequence of 



alternations of states Sq 



si 



. . . s n obtained 

by performing delays di and emitting outputs Oj. Here we 
consider a network of NPTAs with states being of the form 
(£, v). We construct random runs according to Algorithm Q] 
We start from an initial state (£q, vq) and repeatedly concat- 
enate random successor states until we reach the bound c for 
the given observer clock C. Recall that v(C) is the value of 
C in state (£, v), and the rate of C in location £ is R(C)(£). 
We use the notation © to concatenate runs and tail(run) to 
access the last state of a run and delay (fi s ) returns a random 
delay according to the delay density function /j s as described 
in Section HID The statement "pick" means choose uniformly 
among the possible choices. The correctness of Algorithm Q] 
with respect to the stochastic semantics of NPTAs given in 
Section IILTI follows from the Theorem below: 



Theorem 1 Let A be a network of NPTAs. Then: 



f) Statistical Model Checking Algorithms: We briefly 
recap statistical algorithms permitting to answer the following 
two types of questions : (1) Qualitative : Is the probability for 
a given NPTA A to satisfy a property Oc<cf greater or equal 
to a certain threshold ? and (2) Quantitative : What is the 
probability for A to satisfy <>c< c ip- Each run of the system 
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is encoded as a Bernoulli random variable that is true if the 
run satisfies the property and false otherwise. 

g) Qualitative Question.: This problem reduces to test the 
hypothesis H : p = Pa.{ < ^'c<c 1 p) > against K : p < 9. 
To bound the probability of making errors, we use strength 
parameters a and (3 and we test the hypothesis Ho : p > po 
and Hi : p < p\ with po = 9 + Sq and p\ = 9 — Si. The 
interval p — pi defines an indifference region, and p and pi 
are used as thresholds in the algorithm. The parameter a is the 
probability of accepting H when Hi holds (false positives) 
and the parameter (3 is the probability of accepting Hi when 
Ho holds (false negatives). The above test can be solved by 
using Wald's sequential hypothesis testing |[T6l . This test, 
which is presented in Algorithm 12 computes a proportion r 
among those runs that satisfy the property. With probability 1, 
the value of the proportion will eventually cross log(/3/(l — a) 
or log((l — P)/a) and one of the two hypothesis will be 
selected. 

Algorithm 2: Hypothesis testing 

function hypothesis(S:model , ip: property) 

1 r:=0 

2 while true do 
Observe the random variable x corresponding to Oc<c<P 
for a run. 

r :=r + x* log(pi/p ) + (1 — a;) *log((l -pi)/(l -po)) 
if r < log(/3/(l — a)) then accept Ho 
if r > log((l — P)/a) then accept Hi 



end 



h) Quantitative question: This algorithm 1 32 1 computes 
the number N of runs needed in order to produce an ap- 
proximation interval [p — e,p + e] for p = Pr(-0) with a 
confidence 1 — a. The values of e and a are chosen by the 
user and N relies on the Chernoff-Hoeffding bound as shown 
in algorithm [3] 



Algorithm 3: Probability estimation 



function estimate(5':model , ip: property, S: confidence, e. 
approximation) 

1 N := 4*log(l/<5)/e 2 , a := 

2 for i := 1 to N do 
Observe the random variable x corresponding to ip for a 
run. 

a := a + x 
end 
5 return a/N 



V. Beyond "Classical" Statistical 
Model-Checking 

Here, we want to compare pi = f'A.{^c 1 <c 1 fi) and 
P2 = Pa.(Oc 2 <c 2 ( P2) without computing them, with clear 
applications e.g. in determining the possible improvement 
in performance of a new control program. In |fl6l . Wald 
has shown that this problem can be reduced to a sequential 
hypothesis testing one. Our contributions here are (1) to apply 



this algorithm in the formal verification area, (2) to extend the 
original algorithm of lfl6l to handle cases where we observe 
the same outcomes for both experiments, and (3) to implement 
a parametric extension of the algorithm that allows to reuse 
results on several timed bounds. More precisely, instead of 
comparing two probabilities with one common cost bound 
C < c, the new extension does it for all the N bounds i * c/N 
with i = 1 . . . N by reusing existing runs. 

i) Comparison Algorithm.: Let the efficiency of sat- 
isfying Od^ci^i over runs be given by ki — pi/(l — pi) 
and similarly for Oc 2 <c 2 l P2- The relative superiority of "cf2 
over cz>i" is measured by the ratio u = r 2 - — P2< - 1 ~ Pl \ . If 
u = 1 both properties are equally good, if u > 1, <f2 is better, 
otherwise (pi is better. Due to indifference region, we have 
two parameters uo and ui such that uq < ui to make the 
decision. If u < uo we favor (fi and if u > Ui we favor 
(fi2- The parameter a is the probability of rejecting (p-y when 
u < uo and the parameter f3 is the probability of rejecting 
if2 when u > u±. An outcome for the comparison algorithm 
is a pair {xi,x 2 ) = (n (= Cl < Cl ^,r 2 (= Oc 2 <c 2 f2) 
for two independent runs ri and r2- In Wald's version (lines 
10-14 of Algorithm @), the outcomes (0,0) and (1,1) are 
ignored. The algorithm works if it is guaranteed to eventually 
generate different outcomes. We extend the algorithm with 
a qualitative test (lines 5-9 of Algorithm [4]i to handle the 
case when the outcomes are always the same. The hypothesis 
we test is P^(n |= Cl < Cl fi = r 2 \= Oc 2 <c 2 (p2) > 9 
for two independent runs n and T2- Typically we want the 
parameters p' = 9 + So (for the corresponding hypothesis Ho) 
and p[ = 9 — Si (for Hi) to be close to 1. Our version of 
the comparison algorithm is shown in algorithm [4] with the 
following initializations: 



logOi)-log(«o) ' ' 



log(tii)-log(u ) ' 



(ui)-log(ti ) 



Algorithm 4: Comparison of probabilities 



10 

n 

12 
13 
14 



function comprise(S':model , ipi, tfe: properties) 
check := 1, q :— 0, t :— 
while true do 

Observe the random variable Xi corresponding to ip\ for a 
run. 

Observe the random variable X2 corresponding to 4>2 for a 
run. 

if check = 1 then 

X = (ll == X2) 

q = g +:r*log(pi/Po) + (l-a;)*log((l-pi)/(l-p'o)) 
if q < log(/3/(l — a)) then return indifferent 
if r > log((l - P)/a) then check = 

end 

if xi 7^ X2 then 

a = a + c, r — r + c 
if xi = and X2 = 1 then t :=t + l 
if t < a then accept process 2. 
it t > r then reject process 2. 



end 



end 
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j) Parametrised Comparisons: We now generalise the 
comparison algorithm to give answers not only for one cost 
bound c but N cost bounds i * c/N (with i = 1 . . . N). This 
algorithm is of particular interest to generate distribution over 
timed bounds value of the property. The idea is to reuse the 
runs of smaller bounds. When <>c<cfi or <>c<cf2 holds on 
some run we keep track of the corresponding point in cost 
(otherwise the cost value is irrelevant). Every pair or runs gives 
a pair of outcomes (x\,X2) at cost points (01,02). For every 
i = 1 ... N we define the new pair of outcomes (y^ , ) = 
(xi A (i ■ c/N > t\ ■ rate c ),x 2 A (i ■ c/N > t 2 ■ ratec)) 
for which we use our comparison algorithm. We terminate the 
algorithm when a result for every i th bound is known. 

Let a, r, c be the parameters of the previous comparison 
algorithm. Let a',r',c' be the parameters of the qualitative 
check of Section [IV] The procedure is shown in Algorithm [5] 

Algorithm 5: The algorithm for parametrised probabilities 
comparison 

function comprise2(S:model , ipi, ip2'. properties, C: clock, c: 
cost bound, N: # of time intervals) 

1 for i := 1 to N do 

2 I qt :— 0, a'i := a', r\ := r', U := 0, 04 := a, r; := r 
end 

3 repeat 

4 Observe X\ corresponding to tp\ for a run at time t 1 . 
s Observe X2 corresponding to ip2 for a run at time tz, 

6 stop := 1 

7 for i := 1 to N do 

s yi :— xi A i * c/N > ti * ratec 

9 y2 := X2 A i * c/N > t 2 * ratec 

10 if resulti = — 2 then 

11 a'i := a'i + c', r'i — r[ + c 

12 if yi = y 2 then q l ■— q l + 1 

13 if qt < a^ then resulti := 0.5 

14 if qi > r'i then resulti := — 1 
end 

15 if resulti < and y-y ^ y2 then 

16 04 := at + c, n = u + c 

17 if y% = and y 2 = 1 then 

18 I U := U + 1 
end 

19 if ^ < at then resulti '■= 1 

20 if tj > r< then resulti := 
end 

21 if resulti < then stop := 0. 
end 

until stop = 1; 



The results for every i th bound are three-valued: means 
(f2 is rejected, 1 means (p 2 is accepted, and 0.5 means 
indifference. 

VI. Case Studies 

We have extended Uppaal with the algorithms described 
in this paper. The implementation provides access to all the 
powerful features of the tool, including user defined functions 
and types, and use of expressions in guards, invariants, clock- 
rates as well as delay-rates. Also the implementation supports 
branching edges with discrete probabilities (using weights), 



thus supporting probabilistic timed automata (a feature for 
which our stochastic semantics of NPTA may be easily ex- 
tended). Besides these additional features, the case-studies 
reported below (as well as the plots in the previous part of the 
paper) illustrate the nice features of the new plot composing 
GUI of the tool. For more results including models of the 



case-studies see http://www.cs.aau.dk/~adavid/smc/ 

k) Train-Gate Example: We consider the train-gate 
example [ 33 1, where N trains want to cross a one-track bridge. 
We extend the original model by specifying an arrival rate 
for Train i ((i + 1)/N). Trains are then approaching, but 
they can be stopped before some time threshold. When a 
train is stopped, it can start again. Eventually trains cross the 
bridge and go back to their safe state. The template of these 
trains is given in Fig. [5ja). Our model captures the natural 
behavior of arrivals with some exponential rate and random 
delays chosen with uniform distributions in states labelled with 
invariants. The tool is used to estimate the probability that 
Train and Train 5 will cross the bridge in less than 100 
units of time. Given a confidence level of 0.05 the confidence 
intervals returned are [0.541,0.641] and [0.944,1]. The tool 
computes for each time bound T the frequency count of runs 
of length T for which the property holds. Figure f5Jb) shows a 
superposition of both distributions obtained directly with our 
tool that provides a plot composer for this purpose. 



Safe 

(1+id):N*N 




Appr 

a) x <= 20 

0,025 
0,020 
£"0,015 
n 0,010 
g 0,005 



10 

b) 



stop[id] 



ll 



] Train(5) 
| Train(O) 



80 150 
Time 



220 290 



Figure 5: Template of a train (a) and probability density distri- 
butions for Or< t Train(0). Cross and OT<*Train(5). Cross. 

The distribution for Train 5 is the one with higher probab- 
ility at the beginning, which confirms that this train is indeed 
the faster one. An interesting point is to note the valleys 
in the probability densities that correspond to other trains 
conflicting for crossing the bridge. They are particularly visible 
for Train 0. The number of valleys corresponds to the number 
of trains. This is clearly not a trivial distribution (not even uni- 
modal) that we could not have guessed manually even from 
such a simple model. In addition, we use the qualitative check 
to cheaply refine the result to [0.541, 0.59] and [0.97, 1]. 

We then compare the probability for Train to cross when 
all other trains are stopped with the same probability for 
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Figure 6: Comparing trains and 5. 



Train 5. In the first plot (Fig. [6] top), we check the same 
property with 100 different time bounds from 10 to 1000 
in steps of 10 and we plot the number of runs for each 
check. These experiments only check for the specified bound, 
they are not parametrised. In the second plot, we use the 
parametric extension presented in Section [V] with a granularity 
of 10 time units. We configured the thresholds uq and u± to 
differentiate the comparisons atttrj = 1 — e and ui = 1+e with 
e = 0.1, 0.05, 0.01 as shown on the figure. In addition, we use 
a larger time bound to visualise the behaviors after 600 that 
are interesting for our checker. In the first plot of Fig. [6] we 
show for each time bound the average of runs needed by the 
comparison algorithm repeated 30 times for different values 
of e. In the bottom plot, we first superpose the cumulative 
probability for both trains (curves Train and Train 5) that 
we obtain by applying the quantitative algorithm of SectionllVl 
for each time bound in the sampling. Interestingly, before that 
point, train 5 is better and later train is better. Second, we 
compare these probabilities by using the comparison algorithm 
(curves 0.1 0.05 0.01). This algorithm can retrieve 3 values: 
if Train wins, 1 if Train 5 wins and 0.5 otherwise. We 
report for each time bound and each value of e the average of 
these values for 30 executions of the algorithm. 

Table I: Sequential and parallel check comparison. 



e 


0.1 


0.05 


0.01 


sequential 


92s 


182s 


924s 


parallel 


5s 


12s 


92s 



In addition, to evaluate the efficiency of computing all 
results at once to obtain these curves, we measure the ac- 
cumulated time to check all the 100 properties for the first 
plot (sequential check), and the time to obtain all the results 
at once (parallel check). The results are shown in Table U 
The experiments are done on a Pentium D at 2.4GHz and 
consume very little memory. The parallel check is about 10 
times fastej0- In fact it is limited by the highest number of runs 

5 The implementation checks simulations sequentially using a single thread. 



required as shown by the second peak in Fig. [6] The expensive 
part is to generate the runs so reusing them is important. 

Note that at the beginning and at the end, our algorithm 
aborts the comparison of the curves, which is visible as the 
number of runs is sharply cut. 

I) Lightweight Media Access Control Protocol: The 
Lightweight Media Access Control (LMAC) protocol is used 
in sensor networks to schedule communication between nodes. 
This protocol is targeted for distributed self-configuration, 
collision avoidance and energy efficiency. In this study we 
reproduce the improved Uppaal model from 11341 without 
verification optimisations, parametrise with network topology 
(ring and chain), add probabilistic weights (exponential and 
uniform) over discrete delay decisions and examine statistical 
properties which were not possible to check before. Based 
on 051 . our node model consumes 21, 22, 2 and 1 power units 
per time unit when a node is sending, receiving, listening for 
messages or being idle respectively. 



F — I exp-ohain 
Frl exp-ring 
F — I uni-chain 
r A H uni-ring 



27 34 41 
time 

(a) Cumulative probability of collision over time. 

0.90 





F— I exp-chain 
F3 uni-chain 



1.4 2.8 4.2 
collisions 

(b) Probability of having various numbers of collisions. 

Figure 7: Collision probabilities when using exponential and 
uniform weights in chain and ring topologies. 

Fig. [7^ shows that collisions may happen in all cases and 
the probability of collision is higher with exponential decision 
weights than uniform decision weights, but seems independent 
of topology (ring or chain). The probability of collision stays 
stable after 50 time units, despite longer simulations, meaning 
that the network may stay collision free if the first collisions 
are avoided. We also applied the method for parametrised 
probability comparison for the collision probability. The res- 
ults are that up to 14 time units the probabilities are the 
same and later exponential weights have higher collision 
probability than uniform, but the results were inconclusive 
when comparing different topologies. 

The probable collision counts in the chain topology are 
shown in Fig. [7J3, where the case with collisions has a 
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probability of 87.06% and 89.21% when using exponential 
and uniform weights respectively. The maximum number of 
probable collisions is 7 for both weight distributions despite 
very long runs, meaning that the network eventually recovers 
from collisions. 

The probable collision count in the ring topology (not 
shown) yields that there is no upper bound of collision count 
as the collisions add up indefinitely, but there is a fixed 
probability peak at collisions (87.06% and 88.39% using 
uniform and exponential weights resp.) with a short tail up to 
7 collisions (like in Fig. IT]?), long interval of probability and 
then small probability bump (0.35% in total) at large number 
of collisions. Thus chances of perpetual collisions are tiny. 

Fig. [8] shows 
energy consumption 
probability density: 
using uniform and 
exponential weights 
in a chain and a 
ring topologies. Ring 
topology uses more 
power (possibly due 
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Figure 8: Total energy consumption. 

to collisions), and uniform weights 
use slightly less energy than exponential weights in these 
particular topologies. 

m) Duration Probabilistic Automata: Duration Prob- 
abilistic Automata [19] (DPA) are used for modeling job-shop 
problems. A DPA consists of several Simple DPAs (SDPA). An 
SDPA is a processing unit, a clock and a list of tasks to process 
sequentially. Each task has an associated duration interval, 
from which its duration is chosen (uniformly). Resources are 
used to model task races - we allow different resource types 
and different quantities of each type. A fixed priority scheduler 
is used to resolve conflicts. A DPA example is shown in Fig. [9] 

[n=4] [ re = 3] 

• p^l — -C^) — |m 





• K J ' I 2 - 3 ! 

[n = l,re = 2] [n = 2,re = l] 

Figure 9: Rectangles are busy states and circles are for waiting 
when resources are not available. There are r\ = 5 and = 3 
resources available. 

DPA can be encoded in our tool (with a continuous or 
discrete time semantics) or in Prism (discrete semantics), 
see the technical report |36j. In PRISM, integer and boolean 
variables are used to encode the current tasks and resources. 
Prism only supports the discrete time model. In Uppaal, a 
chain of waiting and task locations is created for each SDPA. 
Guards and invariants encode the duration of the task, and an 
array of integers contain the available resources. The scheduler 
is encoded as a separate template. We omit the resources 
and durations from the table for simplicity, they are chosen 
arbitrarily for the experiment. For Uppaal, both a discrete 
and continuous time versions have been implemented. The 
performance of the translations is measured on several case 
studies and shown in Tables [XT] and [HI] In the hypothesis 



testing column, UPPAAL (Upp in the table) uses the sequential 
hypothesis testing introduced in Section [IV] whereas Prism 
uses its own new implementation of the hypothesis testing 
algorithm. In the estimation column, both Uppaal and Prism 
use the quantitative check of Section [IV] but UPPAAL is faster 
due to implementation details. For both tools, the error bounds 
used are a = f3 = 0.05. In the hypothesis test, the indifference 
region size is 0.01, while we have e ~ 0.05 for the quantitative 
approach. The results show that Uppaal is faster than Prism 
even with the discrete encoding, which currently is the only 
fair comparison. 

Table II: Tool performance comparison. 





Parameters 


Estimation 




Hypothesis Testing 


n 


k 


Duration 


Prism 




Vpp c 


Prism 


Upp d Upp c 


10 


10 


4,8 


42.2 


8.7 


6.9 


64.1 


1.0 .3 


10 


10 


8,16 


60.3 


11.3 


7.2 


49.4 


.7 .3 


10 


10 


16,32 


91.8 


13.4 


7.0 


77.1 


.9 .4 


10 


10 


32,64 


126.0 


14.8 


7.0 


65.8 


.9 .3 


10 


10 


64,128 


176.8 


16.3 


7.0 


83.4 


.9 .3 


20 


20 


64,128 




129.4 


52.2 




5.2 1.6 


20 


20 


128,256 




146.4 


52.1 




8.1 1.8 


20 


20 


256,512 




173.8 


52.3 




11.6 1.8 



In the first test, we create a DPA with n SDPAs, k tasks per 
SDPA and no resources. The duration interval of each task is 
changed and the verification time is measured. In the second 
test, we choose n, k and let m be the number of resource types. 
The resource usage and duration intervals are randomised. The 
query for the approximation test is: "What is the probability 
of all SDPAs ending within t time units?". In the verification 
test, we ask the query: "Do all SDPAs end within t time units 
with probability greater than 40%?". The value of t varies for 
each model as it was computed by simulating the system 369 
times and represent the value for which at least 60% of the 
runs reached the final state. 

Table III: Comparison with various durations. 



Parameters 


Estimation 


Hypothesis Testing 


n 


k 


m 


Prism Upp d Upp c 


Prism 


Upp d 


Vpp c 


4 


4 


3 


2.7 1.3 1.0 


2.0 


.1 


.1 


6 


6 


3 


7.7 3.4 2.6 


3.9 


.2 


.3 


8 


8 


3 


26.5 6.9 5.6 


16.4 


.4 


.2 


20 


40 


20 




>300 


34.2 


24.4 


30 


40 


20 




>300 


57.3 


38.0 


40 


40 


20 




>300 


67.4 


70.0 


40 


20 


20 




>300 


40.0 


35.4 


40 


30 


20 




>300 


55.5 


51.4 


40 


55 


40 








219.5 


50 


55 


40 








323.8 


55 


40 


40 








307.0 


55 


50 


40 








342.7 



VII. Conclusion and Future Work 

This paper proposes a natural stochastic semantics for 
networks of priced timed automata. The paper also explains 
how Statistical Model Checking can be applied on the resulting 
model, handling case studies that are beyond the scope of 
existing approaches. 
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The case studies show that models are more expressive, 
the tool is faster and capable of handling larger models than 
the scope of the state-of-the-art model-checker of stochastic 
systems. The extended property language allows quantification 
of events with a limited impact in terms of probability and cost 
complementing critical property checks. Hypothesis testing 
has an order of magnitude advantage in verification time over 
probability estimation, thus provides an opportunity to gain 
leverage when more information is available. 

There are many directions for future research. For example, 
the designer may have some prior knowledge about the prob- 
ability of the property violation. This information could be 
used in a Bayesian fashion to improve the efficiency of the 
test. If the system is assumed to be "well-designed", one can 
postulate that the property under verification should rarely be 
falsified. In this case, the statistical model checking algorithms 
will be efficient to compute the probability of absence of 
errors. Unfortunately, they will not be efficient to compute 
the probability of making an error. We propose to overcome 
this problem by mixing existing SMC approaches with rare- 
event techniques 11371 . Finally, it would also be of interest to 
consider more elaborated properties [38|-[41| or black-box 
systems IH31 . 
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